Effective Storage of Electrons in Water by the Formation of Highly Reduced Polyoxometalate Clusters

Aqueous solutions of polyoxometalates (POMs) have been shown to have potential as high-capacity energy storage materials due to their potential for multi-electron redox processes, yet the mechanism of reduction and practical limits are currently unknown. Herein, we explore the mechanism of multi-electron redox processes that allow the highly reduced POM clusters of the form {MO3}y to absorb y electrons in aqueous solution, focusing mechanistically on the Wells–Dawson structure X6[P2W18O62], which comprises 18 metal centers and can uptake up to 18 electrons reversibly (y = 18) per cluster in aqueous solution when the countercations are lithium. This unconventional redox activity is rationalized by density functional theory, molecular dynamics simulations, UV–vis, electron paramagnetic resonance spectroscopy, and small-angle X-ray scattering spectra. These data point to a new phenomenon showing that cluster protonation and aggregation allow the formation of highly electron-rich meta-stable systems in aqueous solution, which produce H2 when the solution is diluted. Finally, we show that this understanding is transferrable to other salts of [P5W30O110]15– and [P8W48O184]40– anions, which can be charged to 23 and 27 electrons per cluster, respectively.


SI-2: Comparison with different counter cation of W-Dawson cluster
Geometrically Li6 [P2W18O62]. provides 2 kinds of identical W atoms environment in W-Dawson cluster, including 2 belt W6 hexagonal belts in between 2 triangular caps W3O13. In each W6 belt the octahedral are alternatively sharing edges and corners while the 2 belts are joined by 6 common corners. And the electron paramagnetic resonance (EPR) spectroscopy S9 results show a g-factor 1.855 EPR resonance can be observed in the 1e reduced sample which corresponds to the 1e-delocalization electrons on one of the internal hexagonal belts. 1 See main text for further discussion of the EPR spectra of frozen solutions of Li-{P 2 W 18 }-ne (100 mM; n = 1, 2, 3, 4, 5, 6, 12, 17).

EPR fitting:
The EPR spectra were fitted using PHI 3.1.1 program 2 adopting the following effective spin Hamiltonian: The best fits were obtained using the following parameters: g = 1. Na2WO4 solution (δ=0 ppm) according to the IUPAC recommendation. High quality spectra, with good signal-to-noise ratio to observe tungsten satellites, required 1 day of spectrometer time, or until a spectrum of suitable resolution was obtained.  The dc-magnetisation of quartz container (green circles); 0.056 ml of deionised water in the quartz container (red circles) and the magnetisation of 0.056 ml of water without contribution from the container in the temperature range 100 -250 K in a magnetic field of 1000 Oe. The range is consistent with nearly temperature independent magnetism of both water and quartz container.

SI-6: Small-angle X-ray scattering (SAXS) data
Benchtop X-ray scattering data were collected on an Anton Paar SAXSess instrument using Cu-Kα radiation (1.54 Å) and line collimation. The instrument was equipped with a 2-D image plate for data collection in the q = 0.018-2.5 Å -1 range. The lower q resolution is limited by the beam attenuator. Approximately 20 mmolar solutions were measured in 1.5 mm glass capillaries (Hampton Research). Scattering data of neat solvent was collected for background subtraction. Scattering was measured for 30 min for every experiment. We used SAXSQUANT software for data collection and treatment (normalization, primary beam removal, background subtraction, desmearing, and smoothing to remove the extra noise created by the desmearing routine). All analyses and curve-fitting to determine Rg, size, shape and size distribution were carried out utilizing IRENA macros with IgorPro 6.3 (Wavemetrics) software. 1 To simulate scattering data of P2W18 monomer and dimer aggregates, we used SolX software. 2     Table 1 in manuscript for fitting parameters. Red markers is experimental data (plus error bars), gray line is the model fit of the scattering data.  Table 1 in manuscript for fitting parameters. Red markers is experimental data (plus error bars), gray line is the model fit of the scattering data.  Table 1 in manuscript for fitting parameters. Red markers is experimental data (plus error bars), gray line is the model fit of the scattering data.

Electronic structure of Dawson anion and theoretical analysis of initial reduction states
Using a GGA functional (BP86) and Slater TZP basis set, we have computed the redox potentials with respect to NHE. Taking an absolute free energy of -4.24 eV for NHE, 1 the redox potential for the first reduction was found to be -0.04 V, which is about 200 mV more negative than the observed potential for the Dawson anion. In terms of molecular orbital energies and neglecting the bielectronic effects, this would mean that the energy of the LUMO should be somewhat more negative than what we get. In line with this relatively small deviation, when we compare experimental and theoretical redox potentials we reproduce rather well the trends for two successive electron reductions, as shown in Table S8 a) Reduction potentials in V vs SHE; b) reduction potential differences for two successive electron reductions in mV. All the reduction potentials were computed with a BP86 functional and a Slater TZP basis set.
To further prove that reduced {P2W18O62} clusters are prone to associate protons to its framework at highly acidic conditions, we carried out a DFT-MD simulation for the {P 2 W 18 }-4e species. Previous to this, we performed an initial classical MD simulation of {P 2 W 18 }-4e at pH ≈ 1 that showed that hydronium cations tend to accumulate close to the highly negatively charged anion (q = 10-) (see Figure S8-1). Thus, we started the DFT-MD simulation from a configuration in which the POM is surrounded by several H3O + cations that are hydrogenbonded to its terminal oxo groups. Figure S8-2 shows the evolution of the number of protons associated to the cluster during 5.5 ps of simulation. At t = 0.5 ps the system has already taken 3 protons at terminal oxo sites and after more than 2 ps oscillating between geometries with a number of protons ranging from 2 to 5, the system mostly remains tri-protonated until the end of the simulation.  Releasing another proton to give H{P 2 W 18 }-4e is also endothermic although to a lower extent (4.5 kcal mol -1 ). Thus, it might be reasonable to think that at pH close to 1 W 18 -4e is di-or mono-protonated. From this point on and assuming a proton-coupled nature for every following reduction process, it might be sensible to expect a number of protons close to 16 in the structure of the 18 electron-reduced system.   No additional acid/base processes were observed during the remaining ~2 ps of simulation (see      For this reason, we performed partial geometry optimizations to keep only one hydrogen bond in the optimized structure, which is the average number of intramolecular hydrogen bonds obtained from the CPMD simulation with explicit solvent. To do so, we constrained all the P···O-H angles but one to 160⁰ to mimic the geometry that the hydroxo groups adopt in solution to form hydrogen bonds with solvent water molecules instead than with other oxygen atoms of the POM structure. This approach allowed a systematic analysis of different structures with different bridging:terminal ratios, all them with 17 protons (Table S8-2). Comparison of the relative free energies strongly suggests that the hyper−reduced cluster in aqueous solution might combine protons at bridging and terminal oxygen sites, the species with 7 protons at bridging positions and 10 at terminal ones (7:10) being the most likely proton distribution among the analyzed ones (see second column of Table S8-2Table ). Besides, the 4:13 (+3.2 kcal·mol −1 ), the 10:7 (+10.2 kcal·mol −1 ) or the 13:4 (+14.2 kcal·mol −1 ) ones and probably those presenting a proton distribution in between are also rather likely considering that most of the POMs in solution are not as isolated monomers but forming part of a supramolecular agglomerate, in which intermolecular interactions could alter their relative stabilities.

Polyoxometalate aggregation as a function of the reduction state.
Nevertheless, it is important to note that all of them can form agglomerates in a similar manner S35 as shown in Table S8-2 and Figure 4b of the main text and therefore, none of them is incompatible with the proposed mechanism for the super−reduction process. a All geometry optimizations were performed constraining the P···ObridgeH angles to 160⁰ but one, according to the number of intramolecular hydrogen bonds in the CPMD simulation with explicit solvent. The selected phosphorus is that of the POM hemisphere in which the hydroxo group is found. No constrains were used for the 0:17 isomer as it does not contain any proton at a bridging position. Relative Free energies are given in kcal·mol −1 and orbital energies in eV. The integration distances for both RDFs are given in Å and refer to the centers of mass of both species. b Simulation at a POM concentration of 20 mM.
Similarly to what we observed for H{P 2 W 18 }−4e in Figure 3c, the agglomeration of    +18.6 a All geometry optimizations were performed constraining the P···ObridgeH angles to 160⁰ but one, according to the number of intramolecular hydrogen bonds in the DFT-MD simulation with explicit solvent. The selected phosphorus is that of the POM hemisphere in which the hydroxo group is found. The energy of the highest SOMO is given in eV and relative Gibbs free energies in kcal·mol −1 .
Previous electrochemical studies suggested that the six−electron−reduced metatungstate [HnW12O40] (14−n)− exhibits three W-W bonds within a triad that gather all the extra electrons. 4,5 In addition, the formation of the W-W bonds was proposed to be coupled to the transformation of the terminal oxo ligands into more labile aqua ligands. 5 To explore this possibility, we computed two additional structures containing 1 and 6 W-W bonds at the cap regions, labeled as H 17 {P 2 W 18 }−18e (6b:7t:2aqua) and H 17 {P 2 W 18 }−18e (3b:2t:6aqua), respectively. As depicted in Figure S8- Table S8 (Table S8-

Static DFT calculations
To determine redox potentials, we initially optimized the structures at BP86 1 level of theory with a Slater−type TZP basis set using the ADF2016.01 code. 2 Solvent effects of water were taken into account with the COSMO 3 continuum solvent model (ε = 78.4) and relativistic effects were considered by means of the ZORA approximation. 4 The remaining calculations were performed with Gaussian09 (rev D.01) 5 using the hybrid B3LYP functional 6 and a LANL2DZ basis set 7 for all atoms, with the corresponding pseudopotential for W and P centers.
In this case, solvent effects were simulated by means of the IEF−PCM solvent model 8 and frequency calculations were performed at 298 K and 1 atm. To evaluate protonation energies, we used the experimental value of solvation free energy of a proton in solution 9 and the free−energy change in (de)protonation processes was corrected with the standard state correction of 1.9 kcal·mol −1 that accounts for the change in free energy in going from gas phase at 1 atm to the standard state of 1 M in solution at 298 K.

Classical MD simulations.
The systems were simulated by classical MD using the GROMACS 4.5.4 software 10 and the AMBER99 force field. 11 Force field parameters for POMs were obtained following the procedure of Bonet−Avalos, Bo, Poblet et al, 12 using CHELPG atomic charges derived from the electrostatic potential computed in the gas-phase. They were obtained with the Gaussian09 package 5 at the DFT level (BP86 functional) 1 using the LANL2DZ basis set 7 Solvent effects were included in geometry optimizations by using the IEF−PCM model 8 as implemented in Gaussian09. 5 The set of Lennard−Jones parameters for W and O atoms were taken from previous work, 13 Water was represented with the TIP3P model. 14 All simulations were performed with 3D−periodic boundary conditions using an atom cutoff of 14 Å for 1−4 van der Waals and of 10 Å for 1−4 Coulombic interactions and corrected for long−range electrostatics by using the particle−particle mesh Ewald (PME) summation method. 15 The simulations were performed at 300 K starting with random velocities. The temperature was controlled by coupling the system to a thermal bath using the Velocity-rescaling algorithm 16 with a relaxation time of 0.5 ps to keep the NVT canonical conditions throughout the simulation.
Newton equations of motion were integrated using the leap−frog algorithm, 17 and a time step of 1 fs. The bonds with hydrogens were restrained using the LINCS algorithm. 18 Before the production runs, the systems were equilibrated with 5000 steps of energy minimization In all cases, 50 POM anions were embedded in a cubic solvent box of 94.0 3 Å 3 , as well as 300 Li + and the number of H3O + required to neutralize the charge of the system. Force field parameters for hydronium cations were taken from the AMBER99 11 force field and the TIP3P water model. 14 Equilibrium distances and angles were taken from the DFT−optimized geometry, obtained at BP86 1 level with a 6−311G(d,p) basis set 19 for all atoms using Gaussian09 5 and the IEF−PCM 8 continuum solvent model (ε = 78.4). Atomic charges were computed from a single−point calculation in gas−phase 12 using the same level of theory, as done for all the simulated POM structures. This procedure was already followed by Wippf and co−workers to model the agglomeration of polyoxotungstates with H3O + cations as linkers. 20 Simulations at lower POM concentration (20 mM) were performed in a simulation box of the same dimensions but decreasing the number of POMs from 50 to 10 and the number of counter cations, accordingly.

DFT-MD simulations
DFT-MD simulations were using the CPMD 4.1 code, 21 and adopting the generalized gradient−corrected BLYP exchange−correlation functional. 22 The electronic structure was described by expansion of the valence electronic wave functions into a plane−wave basis set, which is limited by an energy cut−off of 70 Ry. The interaction between the valence electrons and the ionic cores was treated using the pseudopotential (PP) approximation.
In the simulation with H 18 {P 2 W 18 }−18e, a number of cations to neutralize the overall charge was included in the periodic box, while in other simulations, the number of cations corresponds to the average coordination numbers within 20 Å from the POM center of mass computed from 5 ns classical MD. Previous classical simulations were performed using 6 Li + because of the 6:1 stoichiometry of the initial POM salt and when simulating low pH, the number of H3O + ions to reach a concentration equivalent to a pH < 1 was added in the box.